###与其他不同之处在于，此代码从TWAS data中筛出来不显著的数据与case和control进行交集，用的是2018sci附件4中的TWAS数据

setwd("E:/0 公共数据库差异情况/db_for_5hmc/TWASdb/")
library(ChIPpeakAnno)
library(meta)
fn=c("SCZ.TWAS","ASD.TWAS","BD.SCZ")

case=read.table("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_p/bias_AShM_BF_no_motif.txt",head=T,sep="\t")
case=case[case$BF_in_DC>10,]
case_t=data.frame(chr=case$chr,start=as.numeric(case[,3]),end=as.numeric(case[,3])+1)
bed_case=toGRanges(case_t,format="BED",header=TRUE)

con.file=c("nobias_AShM.txt","nobias_AShM_total_reads10.txt","nobias_AShM_p.5.txt","nobias_AShM_total_reads10_p.5.txt")
col_names=c("TWAS.db","con.file.name","case_in_region","case_not","con_in_region","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=8))
names(result)=col_names
for(i in 1:length(fn)){
file=read.xlsx("aat8127_Table_S4.xlsx",colNames = T,sheet = fn[i])
file$start=format(as.numeric(file$start),scientific = F)
file$stop=format(as.numeric(file$stop),scientific = F)
file$TWAS.P=as.numeric(file$TWAS.P)

file_sig=file[file$TWAS.P<0.05,]
file_sig=file_sig[!is.na(file_sig$CHR),]
file_sig$unitID=paste0("chr",file_sig$CHR,":",file_sig$start,":",file_sig$stop)
file_sig=file_sig[!duplicated(file_sig$unitID),]

file_t=data.frame(chr=paste0("chr",file_sig$CHR),start=format(as.numeric(file_sig$start),scientific = F),end=format(as.numeric(file_sig$stop),scientific = F))
bed_file=toGRanges(file_t,format="BED",header=TRUE)

file_no_sig=file[file$TWAS.P>0.05,]
file_no_sig=file_no_sig[!is.na(file_no_sig$CHR),]
file_no_sig$unitID=paste0("chr",file_no_sig$CHR,":",file_no_sig$start,":",file_no_sig$stop)
file_no_sig=file_no_sig[!duplicated(file_no_sig$unitID),]

file_t_no_sig=data.frame(chr=paste0("chr",file_no_sig$CHR),start=format(as.numeric(file_no_sig$start),scientific = F),end=format(as.numeric(file_no_sig$stop),scientific = F))
bed_file_no_sig=toGRanges(file_t_no_sig,format="BED",header=TRUE)

ol1=findOverlapsOfPeaks(bed_case,bed_file)
tcase=as.data.frame(ol1$peaklist$`bed_case///bed_file`)
ol1.2=findOverlapsOfPeaks(bed_case,bed_file_no_sig)
tcase2=as.data.frame(ol1.2$peaklist$`bed_case///bed_file_no_sig`)
for(j in 1:4){
con=read.table(con.file[j],head=F,sep="\t")
con_t=data.frame(chr=con[,1],start=as.numeric(con[,2]),end=as.numeric(con[,3])+1)
bed_con=toGRanges(con_t,format="BED",header=TRUE)

ol2=findOverlapsOfPeaks(bed_con,bed_file)
tcon=as.data.frame(ol2$peaklist$`bed_con///bed_file`)
ol2.2=findOverlapsOfPeaks(bed_con,bed_file_no_sig)
tcon2=as.data.frame(ol2.2$peaklist$`bed_con///bed_file_no_sig`)
result_tmp=data.frame(matrix(NA,1,ncol=8))
names(result_tmp)=col_names

a=dim(tcase)[1]
b=dim(tcase2)[1]
c=dim(tcon)[1]
d=dim(tcon2)[1]
result_tmp[,1]=fn[i]
result_tmp[,2]=con.file[j]
result_tmp[,3]=a
result_tmp[,4]=b
result_tmp[,5]=c
result_tmp[,6]=d
result_tmp[,7]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$estimate
result_tmp[,8]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$p.value
result=rbind(result,result_tmp)
}
}
result=result[-1,]
write.csv(result,"TWAS_enrichment_0.01.csv",quote = F,row.names = F)

metaor3<-metabin(case_in_region,case_not,con_in_region,con_not,data=result,sm="OR",studlab = TWAS.db) 
forest(metaor3)